require(rio)
require(tidyverse)
require(cowplot)
require(gridExtra)
require(dbscan)
require(igraph)
require(tidygraph)
require(ggraph)
require(pals)

Loading required package: cowplot


Attaching package: ‘cowplot’


The following object is masked from ‘package:lubridate’:

    stamp


Loading required package: gridExtra


Attaching package: ‘gridExtra’


The following object is masked from ‘package:dplyr’:

    combine


Loading required package: dbscan


Attaching package: ‘dbscan’


The following object is masked from ‘package:stats’:

    as.dendrogram


Loading required package: igraph


Attaching package: ‘igraph’


The following objects are masked from ‘package:stats’:

    decompose, spectrum


The following objects are masked from ‘package:lubridate’:

    %--%, union


The following objects are masked from ‘package:dplyr’:

    as_data_frame, groups, union


The following objects are masked from ‘package:purrr’:

    compose, simplify


The following object is masked from ‘package:tidyr’:

    crossing


The following object is masked from ‘package:tibble’:

    as_data_frame


The following object is masked from ‘package:base’:

    union


Loading required package: tidygraph


Attaching package: ‘tidygraph’


The following object is masked from ‘package:igraph’:

    groups


The following object is masked from ‘package:stats’:

    filter


The following object is masked from ‘package:rio’:

    convert


Loading required package: ggraph

Loading required package: pals

avgs = as_tibble(import("../Genetics/G25/Data/TXT/Global25_PCA_pop_averages_scaled.txt"))
#avgs %>% sample_n(10)
modavgs = as_tibble(import("../Genetics/G25/Data/TXT/Global25_PCA_modern_pop_averages_scaled.txt"))
#modavgs %>% sample_n(10)

str_extract(avgs$V1,"^\\w+?(?=_)") %>% unique(sort = TRUE)

str_extract(modavgs$V1,"^\\w+?$") %>% unique(sort = TRUE)

str_extract(modavgs$V1,"^\\w+?$") %>% unique(sort = TRUE) %>% as.character() %>% paste0(collapse = "|")

avgs = avgs %>% dplyr::filter(!str_detect(V1,"^(ARG|AUS|BHS|BLZ|BRA|BWA|CAN|CHL|CHN|CMR|COG|CUW|DOM|ETH|GUM|HTI|IDN|IND|JPN|KAZ|KEN|KGZ|KOR|LAO|MEX|MNG|MWI|MYS)"))
avgs = avgs %>% dplyr::filter(!str_detect(V1,"^(NPL|PAK|PAN|PER|PYF|Saka|Sarg|SDN|THA|TJK|TKM|TON|TWN|TZA|UGA|USA|UZB|VEN|VNM|VUT|ZAF)"))
avgs = avgs %>% dplyr::filter(!str_detect(V1,"_o$"))
avgs = avgs %>% dplyr::filter(!str_detect(V1,"low"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Abazin|Abkhasian|Adygei|Aeta|Afrikaner|Agta|Ahiska|Akha|Akhvakh|Alawite|Alevi|Algerian|Altaian|Amerindian|Ami|Andian|Arain)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Armenian_[AGHP]|Arora|Assyrian|Asur|Atayal|Australian|Awan|Aymara|Azerbaijani_R|Basque_[ABGLR]|Batak|Bedzan|Bengali|Berber_[AM])"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Arora|Assyrian|Asur|Atayal|Australian|Awan|Aymara|Azerbaijani_R|Batak|Bedzan|Bengali|Bagvalin|Bahun|Baiku|Balti|Baniya|Bantu)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Besermyan|Bhumi|Biaka|Birhor|Bitonga|Blang|Bolivian|Bonan|Bonda|Brahmin|Brahui|Bulala|Bunt|Burmese|Burusho|Buryat|Cachi|Cambodian)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Cameroon|Chama|Chan|Chen|Chipewyan|Chopi|Chukchi|Chuvash|Circassian|Cochin|Colla|Cree|Dai|Damai|Darginian|Datog|Daur|Dharkar|Dinka)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Dolgan|Dong|Dungan|Dusadh|Dusun|Elmolo|Eritrean|Esan_Nigeria|Eskimo|Ethiopian|Even|Evenk|Ezhava|Ezid)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(French_Au|French_B[ei]|French_[COP]|Fulani|Gadaba|Gagauz|Gambian|Ganguela|Gelao|Georgian_[GIJKLMRST]|Gond|Greek_[ADEIKLMNSTW])"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Gujar|Gurung|Hadza|Hakkipikki|Han|Hawaiian|Hazara|Hezhen|Hinukh|Hmong|Ho|Htin_Mal|Hui|Hunzib|Igbo|Igorot|Indonesian|Iranian_[BJLMP])"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Iraq|Irula|Italian_[ABEFJLM]|Itelmen|Jamatia|Japanese|Jarawa|Jat|Jehai|Ju_hoan_North|Juang|Kaba|Kadar|Kaitag|Kalash|Kalmyk|Kamboj)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Irula|Itelmen|Jamatia|Japanese|Jarawa|Jat|Jehai|Ju_hoan_North|Juang|Kaba|Kadar|Kaitag|Kalash|Kalmyk|Kamboj)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Kamma|Kanjar|Karachay|Karaite_Egypt|Karakalpak|Karata|Karen_Sgaw|Karitiana|Kashmiri|Kazakh|Ket|Khakass|Khamnegan|Khanty|Khatri)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Khmer|Kho|Kikuyu|Kinh_Vietnam|Kirghiz|Knanaya|Kohistani|Koinanbe|Kol|Komi|Kongo|Konkani|Korean|Korwa|Koryak|Kosipe)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Kshatriya|Kubachinian|Kumyk|Kurdish|Kurichiya|Kurumba|Kusunda|Kuy_Suay|Lahu|Laka|Lao|Lawa|Lebbo|Lemande|Lezgin|Li|Libyan)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Lithuanian_[RP]|Luhya_Kenya|Luo|Luzon|Macedonian|Mada|Madagascar|Madiga|Magar|Makhuwa|Makrani|Mala|Malay|Manchu|Mandenka|Maniq)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Maniyani_Kerala|Mansi|Manyika|Maonan|Maori|Maratha|Mari|Masai|Mayan|Mbuti|Mende_Sierra_Leone|Miao|Mix|Mlabri|Mogush|Mon|Moroccan_)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Mountain_|Mozabite|Mulam|Murut|Mwani|Nahua|Nair|Nanai|Nasoi|Nasrani|Naxi|Ndau|Negidal|Nenets|Nepali|Newar|Nganassan|Ngumba|Nivkh|Nogai)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(North|Nyah_Kur|Nyan|Ogiek|Onge|Orcadian|Oroqen|Palestinian_|Pallan|Pamiri|Paniya|Papuan|Parsi|Pashtun|Pathan_Bhopal|Piapoco|Pima)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Piramalai_Kallar|Poduval_Kerala|Polish_|Pul|Pumi|Punj|Qiang|Qing|Quechua|Rai|Raj|Ratlub|Reddy|Relli|Rendille|Riang|Rohingya)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Roma_|Romaniote|Ronga|Ror|Rumelia_East|Russian_[BKRSTVY]|Saharawi|Sakha|Sakilli|Salar|Saliya_Kerala|Samaritan|Sandawe|Santhal)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Satnami_Chhattisgarh|Saudi[AB]|Selkup|Sena|Sengwer|She|Shor|Sindhi|Somali|Spanish_[ABCEGLPS]|Sri_Lankan|Sudanese|Surui)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Syed|Tabasaran|Tai_Lue|Tajik|Talysh_Azerbaijan|Tamang|Tamil_Sri_Lanka|Tarkhan_Muslim|Tat|Telugu|Thai|Tharu|Thiyya|Tibetan)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Tikar_South|Tindal|Tlingit|Todzin|Tripuri|Tsez|Tswa|Turkish_[ABDEGKNS]|Turkmen|Tuvinian|Tyagi|Udi|Udmurt|Ulchi|Umbundu)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Uttar|Uygur|Uzbek|Vaniya_Kerala|Velama|Vellalar|Vishwakarma_Kerala|Vizayan|Wa|Welsh|Wichi|Xibo|Yadav_Telugu|Yakut|Yao|Yemenite)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"^(Yi|Yoruba|Yugur|Yukagir|Yukpa|Yuku|Zapotec|Zhuang)"))
modavgs = modavgs %>% dplyr::filter(!str_detect(V1,"_o$"))

avgs$V1 %>% unique(sort = TRUE)
modavgs$V1 %>% unique(sort = TRUE)

myavgs = avgs %>% select(-V1)
row.names(myavgs) = avgs$V1
mymodavgs = modavgs %>% select(-V1)
row.names(mymodavgs) = modavgs$V1

Warning message:
“Setting row names on a tibble is deprecated.”
Warning message:
“Setting row names on a tibble is deprecated.”

allavgs = bind_rows(myavgs,mymodavgs)

mynn = kNN(allavgs,k = nrow(allavgs)-1)
mynn

k-nearest neighbors for 995 objects (k=994).
Available fields: dist, id, k, sort

mynn = rownames_to_column(as.data.frame(mynn$id)) %>% as_tibble() %>% mutate(`0` = row_number(),.after = rowname)

mynn %>% sample_n(10)

mynames = mynn %>% select(rowname) %>% mutate(val = row_number())

mynn = mynn %>% pivot_longer(cols = matches("^\\d"),names_to = "rank",values_to = "val") %>% inner_join(mynames,by = "val") %>% select(c(1,4,2))

colnames(mynn) = c("target","population","distance")
mynn = mynn %>% dplyr::filter(distance %in% 1:5)

mynn %>% sample_n(10)

write_csv(mynn,"../Genetics/G25/Data/CSV/top_10_dist_all.csv")

mynn %>% dplyr::filter(str_detect(target,"^Ukrainian.+[^o]$"))

mynn$distance = as.numeric(scale(as.numeric(mynn$distance)))

myg = graph_from_data_frame(mynn)

myg

IGRAPH b80b3be DN-- 995 4975 -- 
+ attr: name (v/c), distance (e/n)
+ edges from b80b3be (vertex names):
 [1] ARM_Areni_C     ->Levant_Megiddo_MLBA_o1
 [2] ARM_Areni_C     ->Armenian_Syunik       
 [3] ARM_Areni_C     ->ARM_LBA               
 [4] ARM_Areni_C     ->Ostrogothic_Crimea_ACD
 [5] ARM_Areni_C     ->Greek_Cappadocia      
 [6] ARM_LBA         ->ARM_MBA               
 [7] ARM_LBA         ->ARM_Lchashen_MBA      
 [8] ARM_LBA         ->RUS_Alan_MA           
+ ... omitted several edges

myg = as_tbl_graph(myg) %>% mutate(c_auth = centrality_authority(weights = distance))
myg

# A tbl_graph: 995 nodes and 4975 edges
#
# A directed simple graph with 1 component
#
# A tibble: 995 × 2
  name               c_auth
  <chr>               <dbl>
1 ARM_Areni_C      0       
2 ARM_LBA          1.16e- 9
3 ARM_Lchashen_MBA 0       
4 ARM_MBA          0       
5 AUT_IA_La_Tene   7.69e- 6
6 AUT_LBK_N        7.31e-11
# ℹ 989 more rows
#
# A tibble: 4,975 × 3
   from    to distance
  <int> <int>    <dbl>
1     1   502   -1.41 
2     1   825   -0.707
3     1     2    0    
# ℹ 4,972 more rows

myg %>% select(name,c_auth) %>% as.data.frame() %>% arrange(desc(c_auth))

options(repr.plot.width = 26,repr.plot.height = 30,repr.plot.res = 400,repr.plot.bg = "gray20")

#myp = myg %>% ggraph(layout = 'unrooted')#,circular = TRUE) 
myp = myg %>% ggraph(layout = 'stress',circular = TRUE) 

myp + geom_edge_fan(aes(color = distance,alpha = after_stat(index),start_cap = label_rect(node1.name),end_cap = label_rect(node2.name)),strength = 1.25,edge_width = 0.15) + 
        geom_node_label(aes(label = name),label.size = 0,size = 1,alpha = 0.6,fill = "gray15",color = "gray75") + #,repel = TRUE,max.overlaps = 10) + 
        scale_edge_alpha('Edge direction') + 
        scale_edge_color_gradient2("distance",low = "cyan",mid = "snow",high = "orange") + 
        theme_void() + 
        theme(plot.background = element_rect("gray20"),legend.title = element_text(color = "snow"),legend.text = element_text(color = "snow"),legend.position = "bottom") +
        guides(edge_alpha = guide_edge_direction()) + 
        coord_flip()

Warning message:
“Using the `size` aesthetic in this geom was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` in the `default_aes` field and elsewhere instead.”

rowname	0	1	2	3	4	5	6	7	8	⋯	985	986	987	988	989	990	991	992	993	994
<chr>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	⋯	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>
VK2020_GreenlandE_VA	767	765	940	914	870	802	772	771	758	⋯	978	559	656	598	554	586	236	237	841	840
VK2020_SWE_Karda_VA	797	19	756	887	168	23	163	175	305	⋯	978	559	656	598	554	586	236	237	841	840
Iberia_East_IA	332	345	342	849	341	851	848	349	890	⋯	559	656	653	554	598	586	236	237	841	840
Sarmatian_KAZ_Bisoba	661	666	663	659	665	667	310	660	662	⋯	638	598	656	554	586	978	236	237	841	840
RUS_Sidelkino_HG	625	648	582	649	622	650	780	583	749	⋯	571	838	564	656	236	554	978	237	841	840
BGR_EBA	45	517	276	424	30	42	461	46	411	⋯	559	653	656	554	598	586	236	237	841	840
Spanish_Menorca	964	352	967	963	968	965	944	925	415	⋯	559	653	656	554	598	586	236	237	841	840
Swiss_German	971	888	886	857	856	447	408	305	351	⋯	638	559	656	598	554	586	236	237	841	840
HUN_North_Transdanubia	310	661	666	660	663	667	665	659	658	⋯	638	598	656	978	554	586	236	237	841	840
HUN_Avar_Middle-Late_Danube-Tisza	273	268	594	285	263	563	645	272	643	⋯	137	468	811	174	219	512	384	429	841	840

target	population	distance
<chr>	<chr>	<chr>
Turkish_Trabzon	Greek_Cappadocia	4
IRN_Hajji_Firuz_C	Turkish_Trabzon	4
BelgianB	French_Nord	1
Danish	VK2020_DNK_Sealand_VA	3
BelgianC	Swiss_German	4
Swedish	VK2020_SWE_Skara_VA	4
Russian_Pskov	DEU_MA_Krakauer_Berg	3
HUN_Starcevo_N	HUN_LBK_MN	3
Moksha	Ingrian	3
HRV_Pop_RomanP	SVK_LIA	5

target	population	distance
<chr>	<chr>	<chr>
Ukrainian_Chernihiv	Ukrainian_Rivne	1
Ukrainian_Chernihiv	Ukrainian_Dnipro	2
Ukrainian_Chernihiv	Ukrainian_Sumy	3
Ukrainian_Chernihiv	Ukrainian_Zhytomyr	4
Ukrainian_Chernihiv	Russian_Orel	5
Ukrainian_Lviv	Ukrainian_Zakarpattia	1
Ukrainian_Lviv	Polish	2
Ukrainian_Lviv	Ukrainian_Sumy	3
Ukrainian_Lviv	Ukrainian_Zhytomyr	4
Ukrainian_Lviv	Ukrainian_Rivne	5
Ukrainian_Rivne	Ukrainian_Chernihiv	1
Ukrainian_Rivne	Ukrainian_Zhytomyr	2
Ukrainian_Rivne	Ukrainian_Dnipro	3
Ukrainian_Rivne	Polish	4
Ukrainian_Rivne	DEU_MA_Krakauer_Berg	5
Ukrainian_Sumy	Ukrainian_Chernihiv	1
Ukrainian_Sumy	Russian_Orel	2
Ukrainian_Sumy	Ukrainian_Rivne	3
Ukrainian_Sumy	Ukrainian_Zhytomyr	4
Ukrainian_Sumy	Ukrainian_Dnipro	5
Ukrainian_Zakarpattia	Ukrainian_Lviv	1
Ukrainian_Zakarpattia	Croatian	2
Ukrainian_Zakarpattia	Hungarian	3
Ukrainian_Zakarpattia	Slovakian	4
Ukrainian_Zakarpattia	Slovenian	5
Ukrainian_Zhytomyr	Ukrainian_Rivne	1
Ukrainian_Zhytomyr	Ukrainian_Chernihiv	2
Ukrainian_Zhytomyr	Ukrainian_Sumy	3
Ukrainian_Zhytomyr	Polish	4
Ukrainian_Zhytomyr	Ukrainian_Dnipro	5

name	c_auth
<chr>	<dbl>
England_CA_EBA	1.000000000
Bell_Beaker_England	0.173952611
England_LIA	0.116812306
England_MBA	0.109148862
Corded_Ware_CHE	0.089944309
ISL_Viking_Age_Pre_Christian	0.081853731
Bell_Beaker_NLD	0.051372059
Scottish	0.051015268
SWE_Ollsjo_BA	0.049313911
POL_Chlopice_Vesele_Culture	0.046520461
CZE_Unetice_EBA	0.033693051
Bell_Beaker_CZE_late	0.029849571
England_EastYorkshire_LIA	0.029163725
Scotland_LIA	0.028726403
CZE_Unetice_preC	0.027456141
Scotland_LBA	0.025179970
England_MIA_LIA	0.024486887
England_LBA	0.024417841
VK2020_DNK_Sealand_VA	0.023105356
VK2020_SWE_Skara_VA	0.021617028
England_EIA	0.021563486
CZE_Bilina_BA	0.016005905
VK2020_NOR_North_VA	0.014743177
VK2020_IRL_Eyrephort_VA	0.010155855
Bell_Beaker_Bavaria	0.009698243
Scotland_EIA	0.009426029
English	0.008341500
England_Saxon	0.007196527
SWE_LN	0.006523234
VK2020_EST_Saaremaa_EVA	0.006340086
⋮	⋮
Lebanese_Muslim	0
Maltese	0
Norwegian	0
Ossetian	0
Palestinian	0
Portuguese	0
Russian_Pinega	0
Russian_Pinezhsky	0
Russian_Pskov	0
Saami	0
Sardinian	0
Saudi	0
Sephardic_Jew	0
Slovenian	0
Sorb_Niederlausitz	0
Spanish_Murcia	0
Swedish	0
Swiss_French	0
Tarkhan_Sikh/Hindu	0
Tu	0
Tunisian	0
Tunisian_Berber_Matmata	0
Tunisian_Jew	0
Turkish_Rumeli	0
Turkish_Trabzon	0
Ukrainian_Chernihiv	0
Ukrainian_Dnipro	0
Ukrainian_Rivne	0
Ukrainian_Zhytomyr	0
Vepsian	0